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Higher-order corrections to the relativistic perihelion advance 
and the mass of binary pulsars. 

Maurizio M. D'ElisecQ 

Osservatorio S.Elmo 
Via A.Caccavello 22, 80129 Napoli Italy 

We study the general relativistic orbital equation and using a straightforward perturbation method 
and a mathematical device first introduced by d'Alembert, we work out approximate expressions 
of a bound planetary orbit in the form of trigonometrical polynomials and the first three terms of 
the power series development of the perihelion advance. The results are applied to a more precise 
determination of the total mass of the double pulsar J0737-3039. 



I. INTRODUCTION 

The general relativistic orbital equation for a planet re- 
volving around a star is deduced from the Schwarzschild 
line element 



ds 2 = c 2 jdt 2 - 7 ~ 1 dr 2 - r 2 dfl 2 , 
7 = 1- 2r*/r, dtf = d6 2 - sin 2 6> dip 2 



(1) 



where r* = fi/c 2 and /j, = GM are the gravitational 
radius and the standard gravitational parameter of the 
star, respectively. All remaining symbols have their usual 
meaning for this type of problem. According to the 
geodesic hypothesis, the path followed by the planet, con- 
sidered as a test-body not to disturb the metric, can be 
determined using the time-like Lagrangian 



2L = \ — 



1 12 

c ryt - 



7~ 1 r 2 - r*(0 2 + sin 2 9ip 2 ), (2) 



and a variational principle that uses the functional 

%] = / L (q,<i)dT, q = dq/dr, (3) 

where r is the planet's proper time and the function 
q = q(r) collectively denotes the degrees of freedom of 
a possible generic planetary motion. The path actually 
followed is what makes S[q(r)} stationary, i.e. the func- 
tional derivative 



SS[q(T')} 
Sq(T) 



SL[q(r'),q(r')} 
Sq( T ) 



dr = 0, 



(4) 



is zero when computed on the effective motion. From the 
working point of view it is well known that a cancelation, 
in the integrand, of the functional derivative is equivalent 
to that of the Eulcr-Lagrange derivative 



S_ 
Sq 







d_d_ 
dr dq 



(5) 



and so we obtain four second-order differential equations 
(the Euler-Lagrange equations) determining the sought 
time-like geodesic 

dL d dL 

~d drlF <7 = r, 6>, V3, i, q = dq/dr. (6) 



It is worth noting that from the defining equation of the 
proper time, ds/dr = c, we have 2L = c 2 , namely the 
Lagrangian is a constant of the motion. 

We observe that the Schwarzschild metric, and the La- 
grangian, are invariant under the reflection 



(t,r,9,(p) i-4 (t,r,-K-d,(p), 



(7) 



at the hyperplane 8 = 7r/2, so the mirror image of a 
geodesic curve clearly has the same property. In par- 
ticular, if we consider a geodesic which at r = starts 
within the symmetry hyperplane and is tangent to it, 
it must coincide with the transformed geodesic, since 
the initial values of position and velocity determine a 
geodesic uniquely. These considerations are confirmed 
by the analysis of the Euler-Lagrange equation for 9 



' + 



2f; 



v"" 



sin26> = 0, 



(8) 



which admits the solution 9 = 7r/2 satisfying the initial 
conditions #o = 7r /2, 9q = 0. If we reorient the coordinate 
system so that these conditions are met, the motion takes 
place in the equatorial plane, and we can simplify the 
Lagrangian ([2]) assuming sin 9 = 1, 6 = 



2L = c 2 7 t 2 -7" 1 !- 2 -rV. 



(9) 



The coordinates tp and t are both cyclical: they appear 
in the Lagrangian only in dotted form, and this means 
two conservation laws. For the azimuth tp we have 



d_dL 
dr dip 







dL 

dp 



rp 



h = const., 



(10) 



while for t we have 



d dL „ dL , . . . 

r = => — r = c 7* = K — const. . (11) 

dr dt dt y 

The two constants are related to the angular momentum 
and to the energy, respectively. Insertion of the two in- 
tegrals into Eq. ([9]) leads to 



2£ = 7-^-7-^-^ 

c r z 



(12) 



In place of the Euler-Lagrange equation for r, it is easier 
to derive the radial equation from the integral 2L = c 2 , 
obtaining so after multiplication by 7 



h 2 2fi 2/iV 



c- - - 1 = 0. 



(13) 



We need to cast the path equation into a form which 
clearly displays the fact that we arc dealing with a Ke- 
plerian orbit subjected to small relativistic corrections, 
so we find it convenient to eliminate the variable r and 
introduce the angle ip instead. This is possible since r 
does not enter the equation directly, but only via its dif- 
ferential dr, so wc transform the r-derivativc into a tp- 
derivative, according to the identity d/dr = (h/r 2 )d/dip. 
Moreover, the formulas simplify considerably if we re- 
place r by its reciprocal u = 1/r. Denoting by a prime 
differentiation with respect to <p, we have r = h 2 u' 2 and 
Eq. (fl3|) becomes 



2// 



*„ ,3 



1 



W - -£« - 2rV + to c" - -7T = 



I, 2 



h 2 



(14) 



while with a differentiation with respect to ip wc find 



u'{u" + u-j^ -3rV) 



(15) 



A solution of this equation is u' — 0, or u — const., and 
this means a circular orbit. Ruling out this possibility, 
Eq. (fT5)) requires that 



u" + u= -j^ + 3r*u 2 . 



(16) 



Comparing this equation with the classical Binet's for- 
mula for a particle subjected to a central force of magni- 
tude f(u) in polar coordinates 



u + u 



h 2 u 2 



(17) 



it appears that a particle moving in the Schwarzschild 
field will behave as though it were under the influence 
of an effective Newtonian inverse-square force plus an 
additional fourth-power inverse force 



a ir*h 2 \r 

— -| — 



(18) 



in a framework in which proper time is used as indepen- 
dent variable. We remark that for planetary trajectories, 
in Eq. (|16|1 the terms u and n/ti 2 are comparable, while u 
and 3r*u 2 differ by a factor 3r*u. The maximum value of 
this quantity corresponds to the planet which is nearest 
to the star, so, for Mercury, with r* « 1.476 • 10 5 cm, r w 
5.5-10 12 cm, it is 3r*/r ra 1CT 7 . This implies that Eq. ([16]) 
represents in general an oscillator with a constant forcing 
term and a weak quadratic nonlinearity which affects its 
frequency^. 

The orbital equation (|16p is not the complete solution 
to motion problem, which would require the knowledge 



of the function (p(i), being t the coordinate time, which 
is the time measured by an observer at rest at great dis- 
tance from the origin, and therefore the quadrature of 
the equations 



dr 

dip 



1 
hu 2 



dt 
dep 



jhu 2 



(19) 



but this does not concern us here. 

The importance of Eq. (fj"6]) is due to a variety of rea- 
sons. It originates, without any approximation, from an 
exact solution of Einstein's equation. Because of the lack- 
ing of a time variable, it is entirely geometric, so it can 
be employed to deduce the precise shape of the orbit and 
one of the most important post-Newtonian predictions of 
the theory: the advance of the perihelion of an elliptical 
orbit 7 -. The equation, then, is naturally linked to Newto- 
nian dynamics, that in the description of a planetary mo- 
tion provides already an excellent approximation as well 
as a clear connection with observation. Last, it is simple 
when compared with other post-Newtonian equations of 
planetary dynamics, and therefore its mathematical and 
theoretical limits can be clearly defined. We will use it to 
calculate the orbit of a test particle in the gravitational 
field external to a non-spinning spherical mass eventu- 
ally to the order 1/c" for any arbitrary positive integer 
n, and the corresponding formulas of the periastron ad- 
vance to the same order. If the effects predicted fall in 
the range of the observability, mcasurability or indirect 
determination for those physical systems where the equa- 
tion is applicable, then no doubt they should necessarily 
be taken into account in all specific cases. 

A determination of higher-order terms of the perias- 
tron advance of a binary pulsar by using the second 
post-Newtonian (2PN) method have been effected by 
some auctors^^. Their results based on the 2PN theory 
can be applied to the problem of a test particle under 
a strong gravitational field by letting the mass of one 
pulsar theoretically approach zero. Since one can argue 
about the rigor of this and other methods used to han- 
dle the 2PN problem of motion, wc think that a proper 
analysis of the correctness and of the limits of accuracy 
of these approaches should be based on the agreement 
under some convenient limit with the results exposed 
in the present paper. When considering relativistic ef- 
fects, other post-Keplerian phenomena come into play 
in the motion of bodies. For example, one could insert 
the relativistic force (fl8|) as a perturbing radial acceler- 
ation in Gauss equations for the variations of the Kep- 
lerian orbital elements^. In particular, in the expression 
of the time derivative of the mean anomaly, in addition 
to the radial force it appears explicitly the motion of 
the perihelion. We could cite at this regard two recent 
work o 10 i 11 on the secular advance of the mean anomaly in 
binary systems, which could be easily extended to include 
the higher-order effects calculated in this paper. It would 
also be interesting to compare the results so obtained 
with those computed by means of the post-Newtonian 
Lagrangian planetary equations^, and this could be the 



subject of future work. 



II. THE ITERATION METHOD 

In the following we will express h in terms of elliptic el- 
ements of Newtonian approximation, so h — sfWPi where 
p = a(l — e 2 ) is the semi-latus rectum, a is the semi-major 
axis and < e < 1 is the eccentricity. Although we have 
excluded before the circular orbit, the orbital equation 
encompasses also this possibility. To make its structure 
more apparent and to facilitate the calculations, wc cast 
Eq. (|16[) in dimcnsionless form, which represents a well- 
known problem in mathematical physics 



u= 1 



(20) 



where this time u means p/r and e is the pure number 
3r*/p. Thus, when r = a(l ± e), u = 1 =p e. To find 
an approximate solution to this equation, one could use 
regular perturbation theory, writing u in the form of a 
perturbative expansion in powers of e, but it runs in trou- 
ble here, since we wish an expansion that converges for 
all values of the independent variable p. To do this in 
the most convenient way, wc perform a shift in the origin 
w = u — 1, in order to write the equation as that of the 
perturbed harmonic oscillator 



w" + w = e(l + w) 2 . 
When e = the equation becomes 
w'q + w = 0, 
whose general solution is 

woi'f) = Wo(0) cosp + w' (0) simp. 



(21) 



(22) 



(23) 



We shall adopt the initial conditions wo(Q) = e, w' (0) = 
0, which hereafter we shall denote standard, and so 



w o(v) = ecosp. 



(24) 



This choice of the initial conditions leads to a simplifica- 
tion of the algebra. Solution (|2"3j) has frequency one and 
period 27r. The equation of the orbit is then 



w + 1 1 1 + e cos ip 
r 



P 



P 



(25) 



that is an ellipse, where the angle tp = locates the 
position of the perihelion (because there the function wq 
is maximal), and identifies the direction of the apse line, 
the greatest symmetry axis of the orbit in the plane. 

If we try a straightforward iterative perturbation 
scheme to solve Eq. ([21]) attaching the appropriate sub- 
scripts (the iteration numbers) to its sides, for the first- 
order solution w± we obtain the equation 



H'l 



+ w\ = e(l + wq)' 



(26) 



of a forced pendulum where there is a resonant cos^ 
term. If wc do not want that the pendulum oscillates 
with an ever increasing period {w\ must stay small for 
all values of tp), then the external force is not allowed to 
have a Fourier component with the same periodicity as 
the pendulum itself. Note here and in the following that, 
according to the method of undetermined coefficients, the 
particular solution of an equation of the form 



w 



k w = y X n cosnkp, 



is formally expressed by 



k 2 



X n cos nkp 



k 2 (l-n 2 ) 



V 



X n cos nkp 



£i fc9 ( 1 -» a )' 



(27) 



(28) 



whose second term is singular, but can be regularized 
since, applying the L'Hospital's rule, wc have 



lim 



X n cos nkp X\ p sin kp 



►l k 2 {l-n 2 ) 



2k 



(29) 



with p explicitly present outside the argument of the 
trigonometrical function, and therefore growing without 
limits with time t, since tp is a monotonic function of 
t. Obviously this would destroy the stability of the or- 
bit. Therefore this simple perturbation scheme does not 
work for Eq. (|21[) . Things go differently if we rearrange 
Eq. ([HI) in the form 



w" + kfw = e{l + w 2 ), 



l-2e. 



(30) 



The equation is unchanged, but now we have chosen to 
look at the isolated linear term on the right side as part 
of the unperturbed equation, which is 



Wq + k^WQ = 0. 



(31) 



Equation (|31[) now represents an oscillator whose natural 
frequency is fc e , smaller than 1, and its solution obeying 
to the standard initial conditions is 



wq(p) = e cos k e p. 



(32) 



It follows that, with respect to the Newtonian oscilla- 
tor (|22|) . in the relativistic case, to the lowest order, the 
values of r, which trace out an approximated ellipse, do 
not begin to repeat until somewhat after the radius vec- 
tor has made a complete revolution. Hence the orbit may 
be regarded as being an ellipse which is slowly rotating. 
In particular, the angular advance for revolution of the 
apse line is given by 



Auj = 2tt 



1 



1 = 27re + 0(e 2 



(33) 



This is Einstein's perihelion formula, but it represents 
only the first term of a series development in powers of 
e. Our aim is to compute this series up to order e 3 . Al- 
though the first approximation gives a satisfactory degree 



of accuracy for ordinary planetary problems, and the sec- 
ond one can cope with particular astrophysical situations, 
we push a step further the computation because with lit- 
tle extra work we shall exhaust all conceivable theoretical 
needs before the orbital equation fails to represent the 
relativistic motion of bodies in strong-field situations. 

Let us try again the perturbation scheme on Eq. (|30p . 
We have now 



+ fc e V =e(l + tug), 



1 + 2 e ' 



1 , 

e—e cos2k e ip, 



(34) 
(35) 



that we solve with the standard initial conditions. It is 
to be noticed that this time we have not any more the 
resonant term cos k e (p. This equation is of the type 



»'i 



k 2 w\ = A + Bcos2k e t. 



(36) 



The solution is the sum of the general solution of the ho- 
mogeneous part and of a particular integral of the com- 
plete equation. This integral can be found using for- 
mula (HE]). Thus the solution of Eq. (J36j) is 



(37) 



iv\ = E cos k e (p 



tl _ 2 Bk- 2 
Ak e — cos2fc e (p. 



As W\ is of order e as A,B are, we put k~ 2 = 1, while 
the constant E is determined by the standard initial con- 
ditions to be 



and so we find 
u>i = e I 1 



E 



—e 
2 



1 



(38) 



— e 
3 



e—e cos2k e (p, 
6 



cos k e ip 



which we shall write in abridged form 



"'i 



Ai + (e + Ei) cos k e ip + B\ cos2k e ip, 



(39) 



(40) 



where as a notational aid we agree, here and in the fol- 
lowing, that the subscript number i attached to a capital 
letter representing a coefficient wants to emphasize the 
presence of the factor e\ This solution is periodical and 
bounded for all values of tp, and represents the general 
relativistic first-order deviation from the classical ellip- 
tical orbit. In the next iteration, which should give the 
solution correct to order e 2 



w 2 + k 2 w 2 = e(l + w 2 ), 



(41) 



the right-hand side will present again a resonant term, 
since to this order we get 



w 2 + K w 2 = A + H 2 cos k e (p + B cos 2k e ip + C cos 3k e p, 

(42) 

(43) 



where H 2 = ee(2Ai + B±) = e 



; e(12 + 5e 2 ) 
6 ; 



with A, B, C to be written later. To get rid of the reso- 
nant term H 2 cos k e (p on the right side of Eq. (J42]) we will 
use a device which stems naturally from the logical path 
we followed in writing Eq. (|30[) and that was introduced 
for the first time by d'Alembert to control the plague of 
the secular terms present in the integration of the equa- 
tion of the lunar motion, which is of the same type of 
that we are considering herei In order to suppress the 
unwanted cosine term, we add a counter term and write 
Eq. ([301) in the form 



A- 



El 
e 



w = e(l 



w 2 )- 



H 2 



(44) 



and consequently we consider the approximate equation 



w 2 + k 2 2 w 2 = e(l + w\) w\, 



4 



Ho 



(45) 
(46) 



In the right side, the added term suppresses the resonant 
term, since replacing w\ with e cos fc^ ip (the other terms 
of w\ would give terms of order greater than e 2 ), we ob- 
tain just — H 2 cos k^(p. In the left side the coefficient k 2 
will be diminished by the amount H 2 /e, so that we get 



l-2e 



12 + 5e 2 



(47) 



This way we shall obtain an acceptable solution u>2 of 
Eq. ([41]) to order e 2 and, at the same time, the correc- 
tion to the frequency to the same order. The procedure 
can obviously be repeated until we have reached the re- 
quired approximation degree for both solution and fre- 
quency. The essence of this method was rediscovered 
by Lindstedli^ and its practical application was further 
elaborated by Poincarei^. Equation ([45J) has the form 



w 2 + fc j. W2 = A + B cos 2kg ip + C cos 3k<2ip, 

A = - 



(2 + e 2 ) (3e + e 3 ) 2 



B 



C = 



3e 



-e , 



(> 



e , 



(48) 
(49) 

(50) 

(51) 



and its general solution, by Eq. ([28J) , is 



w 2 = E cos fc e 2 ip + Ak 2 ' 



Bk' 



■ COS2fc e 2(y9 



8 



cos 3kg <p, 



(52) 



where now it suffices to put fc^ = 1 + 2e. We thus 
find, by determining E by means of the standard initial 
conditions, the second-order approximation to the orbit 



w 2 = A\ + A 2 + (e + Ex + E 2 )coskgip 

+ {B\ + B 2 )cos2k e 2ip + C2COs3fc, 



;#ip, (53) 



where 



the average advance rate is 



A 2 =e 
E 2 =e 
B 2 = 



2 6 - 3e + 3e 2 - e 3 



29e 3 - 96e 2 + 96e - 288 
144 ' 

2 e 3 -3e 2 + 3e 



C 2 =e 2 ' 



48 



Let us consider now the next equation 



w 3 + k#W3 = e(l + w 2 ) w 2 



(54) 



Since we are interested only in the resonant term of 
0(e 3 ) we do not solve this equation, but we can extract 
quickly from the right side the secular generating term 
H% cos k^ ip (see Appendix) , where 



H 3 = ee(2A 2 + B 2 ) = e 



3 e(36-15e + 15e 2 -5e 3 ) 



(55) 



Following d'Alcmbert's method, this term will be can- 
celed writing Eq. (j5~4")l in the form 



»3 + fc> 3 = e(l + w£)- 



H 2 + -ff 3 ,2 ,2 H 3 
W 2 , K t 3 = k l 2- 



( 6 56) 



and wc have 



fci = 1 - 2e 



12 + 5e 2 2 36- 15e + 15e 2 -5e 3 



(57) 



Further, we find 



kg = 1 — e 



18 + 5e 2 , 126-30e + 45e 2 -10e 



12 



36 



(58) 
, i 30 + 5e 2 2 270-30e + 75e 2 -10e 3 , 

^ 12 36 

(59) 

The rotation for revolution of the periapse, by denoting 
with u its angular measure, is given, to order e 3 , by 

Auj = 2Tr(k^ 1 - 1) =27re+57rn+-e 2 )e 2 

+ 5^3-i e + ^ 2 -i e 3 y. (60) 

This formula agrees with the results obtained with the 
Poincare-Lindstedt method and some of its equivalent 
modifications^^. By denoting with P the anomalistic 
period, that is the time that elapses between two pas- 
sages of the object at its perihelion expressed in days, 



Aw 
~P 



w(rad/d) = H — — m U\ + w 2 + W3 



67T2'* 



a(l-e 2 )P 
157rr* 3 (54- 



157rr* 2 (6 + e 2 ) 
h 2a 2 (l-e 2 ) 2 P 
6e + 15e 2 - 2e 3 ) 



2a 3 (l-e 2 ) 3 P 



(61) 



What is the meaning of e in the solution w n of order 
e"? In the zero-order solution wq, e is the eccentricity, 
but in the successive approximations it loses this charac- 
terization: the symbol e is simply a constant in the open 
interval (0,1) that one introduces in the initial conditions 
to express the shape of the orbit, and that coincides with 
the eccentricity of the osculating Kepler ellipse to the 
path followed by the planet when ip = 0. 



III. APPLICATIONS 

In astrophysical applications Eq. ([BTjl . written in the 
form 



f(r*,a,e,Au) = Q, 



(62) 



is an implicit relation between dynamical and orbital pa- 
rameters characterizing the system under consideration, 
and it can be used to calculate any of them, once known 
the others. Thus, for example, in the solar system r* ,a,e 
are known, and we calculate Aw, the perihelion shift of 
the planetary orbits. It is also evident that a periastron 
advance is highly enhanced by small a's (whence short 
periods) and high orbital eccentricities. A word of cau- 
tion is needed here: we must consider the fact that a very 
large eccentricity can also mean a very small periastron 
distance a(l — e), and so wc must stop precisely at the or- 
bit that just grazes the surface of the star. The problem 
of detecting the motion of periastron (or of apoastron) of 
some highly elliptic extrasolar planets has been consid- 
ered by some authors^^, and of course now the question 
is to determine the magnitude of the higher orders effects 
for some plausible orbital parameters. In the first column 
of Table I for comparison purposes we have inserted the 
data concerning the Sun-Mercury system, while the sec- 
ond and third columns are referred to two hypothetical 
exoplanets of a solar-mass star with great eccentricities 
and/or small radial distances. While is doubtful, given 
the particularly high levels of observational accuracy re- 
quired, that it is actually possible to find planets with or- 
bital characteristics fitted for this purpose, in meantime 
one can imagine a verification of the high-orders peri- 
helion formula achieved by means of a man-made solar 
probe in a carefully planned celestial mechanics experi- 
ment. However, we believe that the higher order terms 
in Eq. (|6T|) must be taken into account in all attempts 
to detect the Sun's Lensc-Thirring effect on the perihelia 
of the inner planets in order to separate with improved 



TABLE I: Mercury and two 


hypothetical 


exoplanets 


System — > 


Sun/Mercury 


Star/ Alpha 


Star/Beta 


M 


1.00 


1.00 


1.00 


a (cm) 


1.475 ■ 10 5 


1.475 ■ 10 5 


1.475 • 10 5 


a (cm) 


5.791 ■ 10 12 


5.791 ■ 10 12 


8.788 ■ 10 10 


e 


0.2056 


0.95 


0.20 


e 


1.038 ■ 10~ 5 


1.038 ■ 10" 5 


1.038 ■ 10~ 5 


P(d) 


87.9 


87.9 


0.164 


tii (rad/d) 


5.703 ■ 10 -9 


5.602 ■ 10 -8 


2.001 • 10" 4 


o;2 (rad/d) 


1.097 • 10" 15 


1.262- 10" 13 


2.652 • 10" 9 


cJ3 (rad/d) 


2.46 ■ 10" 22 


2.873 ■ 10" 19 


4.098 • 10" 14 


Awi (arcsec/yr) 


0.429 


4.220 


1.514 • 10 4 


Ao;2 (arcsec/yr) 


2.26 • 10~ 10 


9.51 ■ 10" 6 


1.998 • 10 _1 


AW3 (arcsec/yr) 


5.09 • 10" 17 


2.16 -10 -11 


3.088 • 10" 6 



we get 



accuracy these two important consequences of general 
relativity. 

Some close binary systems are good subjects for a ver- 
ification and an application of the formulas^ 2 -. Before 
applying Eqs. (|6TJ)) . (j6"Tj) to a concrete example, we re- 
call that a generalization of Enstein's perihelion formula 
says that in a system of two non negligible spherical 
masses, the two-body problem, the gravitational param- 
eter [x means GM = G{m\ + rri2)r&. In this instance 
r* could be loosely named the equivalent or the nominal 
gravitational radius of the system. This result is only 
the first-order approximation to a full relativistic treat- 
ment of the center of mass of two comparable bodies, a 
hitherto unresolved problem, but however we will fully 
exploit the effect of this approximation. For these sys- 
tems, the semimajor axis and the total mass M are not 
directly determinable, while eccentricities, periods and 
periastron advance rate are, so we can use Eq. ([62]) to 
find r*, once we have replaced a with the period P and 
the total mass M (the degeneracy of the mass is removed 
by measurement of another relativistic orbital character- 
istic, the Einstein parameter) using Kepler's third law in 
the form 



i \ 5 / 3 / 
2tt\ It 



2/3 



T 7 f(e) + (p 



2n 



7/3 



4/3 

7) 9{e) 



where 



/(e) = 

9(e) 

He) 



He), 



(64) 



3 


(1-e 2 )' 
15(6 + e 2 ) 



2(1 -e 2 ) 2 ' 
15(54 - 6e 



15e 2 - 2e 3 ) 



4(1 -e 2 ) 



2YA 



To obtain directly the value of the mass M in units of 
the solar mass M@, we can replace in Eq. (|64|) the ratio 
r*/c with the product T Q M, where T ee GM q /c s is the 
mass of the Sun expressed in units of time. In the cur- 
rent literature the expression of Co lacks of the last two 
term of Eq. ([64]) .— Inserting the values of e, P, u deter- 
mined for a given binary system, Eqs. (|64p, (jfiip can be 
numerically solved for r* and M respectively, and will 
give the value of the gravitational radius of the system 
or the total mass. Once known r*, we can complete the 
calculation and use Eq. (|6"Tj) to compute the value of a. 
The value of Co that one determines is due, in general, 
to relativity plus extra classical terms, as the gravita- 
tional quadrupole moment induced by rotation and tidal 
deformations. But strongly self-gravitating objects as bi- 
nary pulsars have a mass pointlike behavior, and thus the 
motion of their periastron must be entirely ascribed to 
relativity. Strictly speaking, since this bodies are rapidly 
spinning, one should consider also the Lense-Thirring ef- 
fect or use the Kerr solution to the Einstein equation, 
but this is an argument for a further study. 

We apply the theory to the double pulsar J0737-3039. 
Its orbital period is the smallest so far known for such 
an object, and it can be determined with great precision 
along with Co. For a such system cumulative effects add 
rapidly, and this allows a meaningful application of the 
formulas, despite the rather small eccentricity. In Tabic 



'P 2 GAf\ 1/3 



Pc 

2^ 



2/3 



*l/3 



(63) 



To be precise, this law is not exactly true in a regime 
of strong fields and little distances, but it can be used 
harmlessly in a perturbativc calculation for the systems 
we are considering within the same constraints we have 
set for the center of mass, that is a first-order approxi- 
mation. Incidentally, this will free our results from other 
complex and still poor understood effects and, besides, 
this is in line with the regime of validity of the equation 
of motion which we have supposed. From Eqs. (|6ip . (|63p 



TABLE II: Pulsar J0737-3039 



e 0.0877775 

P 0.10225156248(d) 

6j 16.89947° (yr" 1 ) <■ 



- 8834.534991(s) 
9.346445651 ■ 10~ 9 (s~ 1 ) 



c 2.99792458- 10 10 (cm/s) 
T 4.925490947 • 10~ 6 (s) 
yr 3.15576 ■ 10 7 (s) 



II we have indicated the determined^ 2 - values of e,P,ui 
and, for the reader's convenience, some numerical values 
used in the calculations, while in Table III are indicated 
all parameters that can be derived through an application 



TABLE III: Pulsar J0737-3039 derived parameters 



To order: 


e 


e 3 


r* (cm) 


3.82014 • 10 5 


3.8199525 • 10 5 


M(M ) 


2.587075 


2.586948 


a (cm) 


8.788391 • 10 10 


8.788680 ■ 10 10 


e 


1.314166- 10~ 5 


1.314057 • 10" 5 


Wl 


16.89947°yr~ 1 


16.89891408°yr" 1 


W2 


- 


0.00055589V- 1 


tJ3 


- 


0.00000002°yr _1 


CJ 


16.89947° yr- 1 


16.89946999 V -1 



of our formulas. In the first column, the computation are 
done to order e, assuming the validity of general relativ- 
ity thorough the use of Einstein's precession formula wi 
of Eq. (j6"Tj) , while in the second column are inserted the 
values we have computed to order e 3 , in particular how 
much of u> comes from Wj, i = 1, 2, 3. Since from the ob- 
servational point of view the single higher-orders contri- 
butions to the precession rate are inextricably combined, 
and so hidden to a direct measurement, we can obtain 
an indirect verification turning to the main consequence 
of this approximation: a diminution by a small amount 
of the total mass of the system with respect to the cur- 
rently accepted value, with the consequent redefinition of 
other parameters, in particular of a, and so we can re- 
fine the model of the system with one that fits the found 
variations. 

Ever more accurate determinations of ui might enable, 
for this as for any other similar system, to test the use- 
fulness of the third-order approximation to the pcrias- 
tron secular motion deduced from the equation of mo- 
tion within the limits of approximation that we have in- 
troduced. However, the results found can be considered 
strictly true for the motion of a test body in the gravi- 
tational field of a central body of mass equivalent to the 
total mass of the binary system. 



IV. APPENDIX 

By simple considerations of powers and arguments, 
and with the aid of the subscript notation, we can quickly 
find the coefficient of the secular term H n once known 
H%, •• • , -Hn-i and the orbit to order e™ -1 . Here is a 
sketch of how we unearth from the right side of Eq. I[53|) 



e(l + tu§) -w 2 = e(l 

e 



the resonant term H% cos kg ip of order e 3 . 
We consider first the multinomial expression 

ew\ = e[Ai + A 2 + (e + E 1 + E 2 )coskgip 

+ [Bi + B 2 ) cos2kg(p + C 2 cos3kgip} 2 , 

(66) 

together with the following algebraic and trigonometrical 
identities 



e(o + b + c + . . . ) 2 = ea 2 + eb 2 + ec 2 + . 

h 2eab + 2eac ■ 

h 2ebc + . . . , 



(cosnkgip) = -cos2nkgip ■ 



1 

2' 

(ib 



acos(n + Vjkgp ■ b cos nkg(p = — cos k&ip 



(67) 
(68) 

(69) 



Here's the argument: in expanding Eq. (|66|) according 
to formula (|67p we can omit to explicitly writing the 
cosines, since we know that each e, E multiplies coskgtp, 
so as each B and C multiplies cos2fc ( 2(/? and cosik^Lp re- 
spectively. This way we can proceed rapidly by inspect- 
ing more concise expressions. Let us filter now Eq. (|66p 
through the sieve represented by the constraints we have 
imposed to isolate the resonant term ~ e 3 cos k^tp. 

We observe first that the squared terms in Eq. (|57|) 
can be deleted because of Eq. (|68| for n = 1, 2, 3. Next, 
we drop all double products of Eq. (f67j) in which the 
sum of the subscript indices is different from 2 and so, 
after multiplication by e, will survive only the terms of 
order e 3 . Last, when we meet products of two cosines, 
we consider only those in which the arguments differ by 
one, and apply to them Eq. (|6T)|) . At the end we obtain 
the following sum of coefficients of cos kg <p 



2eeA 2 + eeB 2 + 2eAiEi + eE x B 



1- 



(70) 



but the last two terms are erased by the e 3 -cocfficicnts 
arising from the rightmost expression of Eq. (|65| . which 
are the constants — 2eAi and —eB\ times the term 
Ei cos kg <p of w 2 , and thus we finally get 



H 3 = ee{2A 2 + B 2 ). 



(71) 



™f) 



e(2Ai + Bi)w 2 , (65) 



* Electronic address: s.elmo@mail.com 

1 Bertotti B., Farinella P., Vokrouhlicky D., 2003. Physics 
of the Solar System. (Dordrecht: Kluwer). p. 313. 

2 Burgay, M.; D'Amico, N.; Possenti, A.; Manchester, R. N.; 



Lyne, A. G.; Joshi, B. C; McLaughlin, M. A.; Kramer, 
M.; Sarkissian, J. M.; Camilo, F.; Kalogera, V.; Kim, C; 
Lorimer, D. R. An increased estimate of the merger rate 
of double neutron stars from observations of a highly rela- 



tivistic sj/stem.Nature, Volume 426, Issue 6966, pp. 531-533 
(2003). 

3 Calura, M.; Fortini, P.; Montanari, E. Post-Newtonian La- 
grangian planetary equations, Physical Review D (Parti- 
cles, Fields, Gravitation, and Cosmology), Volume 56, Is- 
sue 8, 15 October 1997, pp.4782-4788 

4 d'Alembert, A. Recherches sur Differents Points du Sys- 
teme du Monde. (1754-56) Livre 1 Chap. VI, n.27. 

5 Damour, T., Schaefer, G. Higher-order relativistic perias- 
tron advances and binary pulsars, Nuovo Cimento B, Vol- 
ume 101, No. 2. 

6 D'Eliseo, M.M. 2007 The first-order orbital equation, Am. 
J. Phys. 75, 352 

7 d'Inverno, A. Introducing Einstein's Relativity, Oxford 
Univerity Press, Oxford, 1992. 

8 Do-Nhat T. Full asymptotic expansion of the relativistic 
orbit of a test particle under the exact Schwarzschild met- 
ric. Physics Letters A 238 ( 1998) 328-336. 

9 Iorio, L. Are we far from testing general relativity with the 
transiting extrasolar planet HD 209458b 'Osiris'? New 
Astronomy, Volume 11, Issue 7, p. 490-494, 2006. 

Iorio, L. The post-Newtonian mean anomaly advance as 
further post-Keplerian parameter in pulsar binary systems, 
Astrophys. Space Sci., Volume 312, Numbers 3-4. Decem- 
ber, 2007, pp. 331-335. 
11 Lin-Sen Li, Post-Newtonian effect on the variation of time 



of periastron passage of binary stars in three gravitational 
theories, Astrophys. Space Sci., doi:10.1007/sl0509-010- 
0267-4 

12 Kramer, M.; Stairs, I. H.; Manchester, R. N.; McLaughlin, 
M. A.; Lyne, A. C; Ferdman, R. D.; Burgay, M.; Lorimer, 
D. R.; Possenti, A.; D'Amico, N.; Sarkissian, J. M.; Hobbs, 
G. B.; Reynolds, J. E.; Freire, P. C. C; Camilo, F. Tests of 
General Relativity from Timing the Double Pulsar, Science, 
Volume 314, Issue 5796, pp. 97-102 (2006). 

13 Landau, L. - Lifchitz, E. Theorie du Champ, Editions MIR 
Moscou, 1966, p.415. 

14 Lindstedt, A. Astr.Nachr. Nr 2482, 1883. 

15 Mason, D.P. - Wright, C.J., An Improved Singular Per- 
turbation Solution for Bound Geodesic Orbits in the 
Schwarzschild Metric. General Relativity and Gravitation, 
17ol. 16, No. 2, 1984. 

16 Nayfeh, A.H. Perturbation Methods, Interscience, New 
York, 1973. 

17 Pal, A.; Kocsis, B., Periastron precession measurements in 
transiting extrasolar planetary systems at the level of gen- 
eral relativity, Monthly Notices of the Royal Astronomical 
Society, Volume 389, Issue 1, pp. 191-198. 2008. 

18 Schafer, C, Wex, N. Phys. Lett. A 174 (1993) 196; 177 
(1993) 461 (E). 



